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ABSTRACT 


This thesis presents a comparative study of two 
procedures for generating semi-discrete solutions to the pair 


of simultaneous partial differential equations 
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where C and T are functions of x and 6. £4These methods 


are based upon use of Taylor Series and Chebyshev Polynomials 
respectively. In both cases it is shown that the solutions so 
obtained are weakly stable -- ie., convergent but numerically 


unstable. 


These partial differential equations arise in 
conjunction with an extended numerical treatment of the 
penetration theory model for the diffusion of a gas through a 
non-volatile liquid. The procedures investigated result in a 
solution interval significantly longer than that obtained by 


previous researchers, 
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PREFACE 


From the earliest attempts to predict theoretically the 
effect of a simultaneous chemical reaction upon the rate of gas 
absorption by Hatta [14], idealised models of the gas absorption 
process have been used extensively. This thesis employs one of 
the most commonly used models, namely, the penetration theory 
model. A pair of simultaneous partial differential equations which 
give concentration of gas and temperature at a known depth at a 
known time are derived. Two approaches to the design of computer 
algorithms for producing numerical solutions to these differential 
equations are studied in detail. The approaches are based 
respectively upon Taylor Series and Chebyshev Polynomials. The 
Taylor Series approach, although more easily handled on a formal 
basis, is found to possess certain inherent computational 
disadvantages; in this respect the solution in Chebyshev 
Polynomials is more efficient. For example, the former is far 
more prone to numerical error accumulating through numerous 
multiplication and division operations which are circumvented in 
the Chebyshev Polynomials case through the use of special matrix 
relations. Furthermore, "economisation'' technique is not applicable 


to the Taylor Series approach. 


Both procedures are demonstrated to be weakly stable which 


seriously limits the distance interval upon which the solutions 
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obtained remain acceptably accurate. However, both methods 

increase the length of the time interval for accurate solution 
materially beyond that of previous methods. Further extension 
of these intervals appears contingent upon more sophisticated 
approximations for the operator 9/dt, with an attendant gross 


increase in mathematical complexity. 
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CHAPTER I 


SOLUTION IN TAYLOR SERIES 


1.1. Nature of the Problem 


In this chapter, a set of simultaneous partial differential 
equations for the study of concentration of gas and temperature have 
been derived. The theoretical model chosen for the absorption 
system is deliberately simplified, and most suitable boundary 
conditions are adopted. The system consists of a pure gas in 
contact with a stagnant semi-infinite liquid phase. System 
parameters, such as the kinetics of the reaction, the interfacial 
area for mass transfer, the liquid density, viscosity and 
diffusivity, which are known to vary, have deliberately been 
taken as constants. The liquid is assumed to be non-volatile 
and hence the chemical reaction and temperature are assumed to 


vary continuously with regard to both time and distance . 


1.2. Equations for Concentration and Temperature 


Let an x-axis be chosen so that the absorbing fluid 
occupies the region x > 0. If C(x,@) is the concentration of 
absorbed gas per unit volume at a depth x and time © , then the 


rate of flow of gas across unit area perpendicular to the x-axis 
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is -p 2 where D is the mass diffusivity. The rate of decrease 


of C caused by the chemical reaction is K(T)C where K(T) is 
the reaction-rate coefficient which is a function of temperature 
T which in turm is a function of x and @. The mass transfer 


equation for the gas is, therefore, 


Gree) =D >-K(T)C. 


The concentration of gas is initially zero throughout the absorbing 
fluid. The concentration at the boundary at a subsequent time 6 
will be denoted by £(6). Thus the concentration C must satisfy 


the following boundary conditions : 


(1. 2a) c(0,6) = £(@) , 
(1.2b) C(x,0) = 0, 
(1.2c) C(~,0) = 0. 


The chemical reaction which decreases Ge sat wtheyrate 
K(T)C absorbs heat at the rate AK(T)C, where A is measured in 
BIU per lb. mole. However, any temperature change involves a 
rate of increase of heat of the amount cee where o ds the 


molal density and os is the molal specific heat. If the thermal 


conductivity is k, the heat transfer equation is 
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which may be written in the form 


2 
T 
(1.4) ae - 2242 4Aouime =o, 
ox iL NE aeo 
where a, = k/pC, is the thermal diffusivity. 


a 


lis Ts denotes the constant temperature throughout the 
liquid at time 6 = 0 and g(6@) the temperature at the boundary 
at subsequent time @ , then the boundary conditions satisfied by 


T are as follows: 


(1. 5a) T(0,8) = g(8) , 
(1.5D) T(x,0) = To > 
(1.5¢) T(#,6) = T) 


Equations (1.1) and (1.4) form a pair of simultaneous 
partial differential equations for determination of C(x,@) and 
T(x,8). In the special instance that the reaction-rate coefficient 
is a constant Ko> the single equation (1.1) serves to determine 
C(x,0). Substitution of the solution C(x,6) into the equation 


(1.4) provides a single equation for the determination of the 
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temperature T(x,6). 


1.3. Danckwerts Solution for Constant Reaction Rate and Constant 


Boundary Conditions 


Danckwerts[5] solved equation (1.1) with K(T) = Ky and 
with constant £(6) = Co. His solution involves the complement of 


the error function and is as follows: 


1 1 
C= $C, exp[-(K,/D) *x] exfc[—- -(K,9) 7] 
2(D8) ” 
(1.6) 
+4 Ce exp[(K, /D) 35] ertc| = = +(K8)" ] 
2(De)? 

For unit area of the surface at x = 0, the rate of 
absorption of the gas is -p 2 and it follows from (1.6) 
that this rate is 

Tay) C (DK 4 fc(K 0) 2+ ( K 9)? (ree) 
Gi) At » erfc(K. 7K exp : 


In the absence of any chemical reaction, the term 


iC 
K(T)C is omitted from equation (igi). andethesrates sD 2 


is equal to 


(1.8) ¢ (D/ne) ® ; 
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The ratio of the expressions (1.7) and (1.8) is 


equal to 
4 1 
(Gach 9,08) = (1K 6) ter (K 8) ° + exp(-K 8) 


At time 9 = 0 the value of o, 69) is unity. At any subsequent 
time 6 , the value of o, (8) indicates the factor by which the 
rate of absorption of gas has increased during this time interval 


due to the chemical reaction. 


1.4. Effect of Linear Boundary Condition 


The boundary condition of constant £(6) is clearly 
unrealistic if the overall temperature of the system increases as 
the reaction proceeds. Sullivan [22] has treated boundary 


conditions (1.2a) and (1.5a) of the form 


CiaO8) COs) a, + a,T(0,6) ‘ 


ii 


(1.10b) T(0,86) Ty 2 EC A 
where a,a, a5 and To are constants. However, he has taken 
K(T) ._ to be a constant Ko: Under these circumstances, he 


derived the following formulas for C(x,9) and T(x,8) : 
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The equations (1.11) were integrated numerically by 
Sullivan [22] who then computed 9C/ax at x = 0 by the following 


approximation where x) is a small number. 





-3C(0,8)+4C(x, ,8)-C(2x, 58) 


oax* X=O 2x4 


Values of the ratio >, (8) which corresponds to o, (8) under the 


changed conditions, were found from 


-D(9C/8x) 24 
(1.14) eat epee LAT 
C(0,6) (D/78) * 


and were plotted as functions of M= (nk) 4. The values of 
parameters were chosen to correspond to several physical situations 
in which co, at a pressure of one atmosphere is absorbed into 

1M NaOH. 


The ratio 4 (8) was derived under the assumption of 
a constant reaction-rate K(T) = Ko and a constant surface 
concentration C(0,9) = Co. The ratio >» (8) was derived under 
the assumption that K(T) = Ko and that the surface concentration 
and temperature are both linear functions of time. The study of 
C(x,6@) and T(x,@) when K(T) is not constant is the subject of 


the following sections and the subsequent chapters. 
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1.5. Characteristics of the Two Partial Differential Equations 


To a large extent, the strategy to be adopted for the 
solution of a partial differential equation depends upon the 
nature of the equation, although the initial and boundary 
conditions play very prominent roles. In this section, a study 
of the characteristics of the equations (1.1) and (1.4) is 


made and domains of solution obtained. 


Consider the diffusion equation (1.1), namely 
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ox 


where C, 9C/ax and 9C/36 are given on a curve m in the 


x - 6 plane. Since -C = C(x,6), we may write 
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The equations (1.15), (1.16b) and (1.16c) lead to 
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Hence the characteristic equations in the x - 9 plane, beneath 


which the solution is valid, are determined by 
2 
GIESAn SD Did) sa 00. 
which gives two coincident characteristics, namely, the lines 
§ = constant, 
and hence the equation is parabolic. These lines run from 
x= 0 to x=. Hence the required region of the x - 6 plane 


is enclosed by characteristics and, therefore, our solutions are 


valid. 
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Similar remarks apply to the temperature equation (1.4), 


namely, 
C1520) Cie ey eee gee UE 
e 
p 
and, therefore, we can proceed with the solution of both the equations. 


1.6.. Treatment of Linear Reaction Rate and General Boundary Conditions 


The function “K(T)’ islLa,continuous function of. IT. Provided 


the temperature changes are small, it may be approximated by putting 
(a2 4) KC) 2292s. 


where & and m are suitably chosen constants. It will hereafter 
be supposed that the temperature T is measured as the excess of 
temperature above the constant temperature at time @ = 0. Thus 


T = 0 and equations (1.1) and (1.4) may be expressed in the 


form 
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L425) L=2/D, M=m/D, o = (1/a,) (A/pC,)D = AD/k . 


With K(T) chosen as in (1.21), the solutions C(x,@) 
and T(x,0) of equations (1.22) are also functions of M and 


they may be expanded as infinite series of the type 
(1. 24a) C(x,6) = JC (eM , 
(1.24b) CA) a CTE) 


where the summations are in s from zero to infinity. 


Substitution of the expansions (1.24) into (1.22), 
and comparison of the terms that involve equal powers of M, leads 
to successive pairs of equations for the determination of the 
functions C (x,8) and T (x58) which in turn are successive 
terms in approximations to C(x,@) and T(x,6) respectively. 


For example, the terms independent of M lead to the equations 
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For the linear boundary conditions (1.10), the 
expressions (1.11) are solutions of equations (1.25) but 
they are not suitable for substitution in the differential equation 
(1.26). In the following section, the solution of equations 


(1.25) is obtained in a more convenient form. 


1.7. Semi-Discrete Solution for C(x, 8) and T (x5 8) 
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In a manner similar to that of section 1.5 it can be 
seen that the partial differential equations (1.25), (1.26) and 
(1.27) are of the parabolic type with characteristics which 
satisfy the equation @ = constant. Hence a valid solution may be 


obtained for x > 0, 6@> 0. We therefore rewrite the above 
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Let fy and g. denote £(6.) and g(6.). The 


Potnacyconditions) fox ove and Ti are then 
(1.32a) co) (xe0) = fy 1") (xe0) = g, (¢ > 0) 
(1.32b) Sp Yreaycrepeiger p69) auth 

(1,32c) CO) (xen) = 0. Teme) = 0. 
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For general values of r, the solution of (1.30a) has 
the form 
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and if this expression is substituted into (1.30a), there result 
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(1.34b) B =f 


These recurrence relations may be used to provide the 
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Equations (1.34a) determine BE n 2t points of a net in 
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values of B. 5 in reverse order, starting from By =f 
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which the points of co-ordinates (rt+l,r) and (r,0) form 
boundaries on which the values of B. n are known by virtue of 
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(1.34b) and a similar previous use of (1.34a). The Spd at 
internal nodes of the net may be determined by steps along lines 
of constant. r values. In contrast, if equation (1.25a) were 
integrated directly by the use of a difference net Liebe x, 
plane, then the boundary conditions on x involve two boundaries 
x =0 and x =o. However, in the technique adopted by us, the 
inclusion of the exponential factor in expression (1.33) ensures 


that the boundary condition at x = ~ is satisfied and it results 


in a simplified difference net for the B. a 


For the special case in which K(T) = Ko C(x, 8) 
coincides with. C(x,0) and the ratio $¢(6) analogous to (1.9) 
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It may be remarked in passing that the recurrence 
relations (1.34) are significantly more convenient for numerical 


computation than the integral expressions (1oi2)..) [he recurrence 
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formulas are well-behaved in the sense that they may be repeated 
several thousand times without introduction of serious round-off 


errors. 


The co-efficients B ef given by the relations (1.34) 
S) 
determine the function of) (x) and hence the concentration 


C (x58) at times 6 = rA@ at depth x. 


In a similar manner, the temperature T (x58) may be 


obtained from equation (1.30b). For r = 0, it has for solution 
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For general values of r, the solution 1°) (x) may be 
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Substitution of the expressions (1.33) and (1.37) into 


the equation (1.30b), and the use of relation CG. me = 1, leads 


to the equation 
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and to the following recurrence relations involving the rr : 
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except those of the type Ser . The substitution also leads to 


the equation 


(1.42) Mig oe 


and to the following recurrence relations involving the Sea 
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fors n=O M)r =1, “with Hg =:0. for, S > %. 


These recurrence relations determine all H e in the reverse 
’ 


order in terms of the B © starting with the known Ho =]. 


> 


To sum up, for each value of r > 0, the relations (1.34) 


may be used to determine the Bey oe equations (1.43) may be 
> 
used to determine the M41 nw and the equations (1.41) may be 
> 
used to determine the Ge aitn except for beer a: Equations 


(1.40), (1.42) and (1.39) then serve to determine Lig? 


M and Gy respectively. Initial values corresponding 


reall aC 


t6. i=s.0 ~are 





(1.44a) Bio = f : 
(1.44b) Gio mn lees 
(1.44e) Het ah 
rf, 
(1.44d) Ly = Bi + ey ‘ 
“Af, 
(1.44e) M. 
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1.8. Semi-Discrete Solution for C, (x, 8) 


The semi-discrete method may be used to solve the equations 
ooh) De?) 
(1-26). for C, (x58) = Cc, (x) and T, (x58) = T; (x) 
Representation of the derivatives as in relations (1.29) and 
substitution into the equation (1.26a) leads to the following 


ordinary differential equation for the determination of sey 


g2c (Ft) 

(1.45) 1 a ac sttt) > gc?) As (etl), (ett) 
2 al il fe) fo) 
dx 
In view of the known expressions for ee and eee 
the above equation may be written as 
2_(rt+1) 
cc are 

af Cre ec) 2 n 
(1.46) 5 ~ aC, = BC, + Le] exp (-vx) ) P ti .n® 

dx n=o 

s, 2% “ 
+ M +1 exp (-2a “x) ) ay ee r 
n=o 

where and v are the known constants determined 


Erie: eerien 


by the equations (in which O.< n < 2r - 2) 
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(1.47b) Ce. ) Hehe 
S=0 
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Since the functions of) (x) were chosen to satisfy the 
boundary conditions (1.32), the corresponding conditions on 


Gat) are zero conditions. 
When r = 0, the solution oy” (x) of equation (1.46) 


is as -follows: 


MQ 
30, 


L,P . 
(1.48) cokes = a [exp (-vx) -exp (-a x) ] -- 


y i 
7 [exp (-20 °x)-exp (-a *x) ] 
v=o 





For general values of r, the solution may be expressed 


in the form 


r-l 2r-2 


(xr) a bs n ~ n 
(1.49) Cy (x) = exp(-a “x) ) Un* + exp (-vx) ) Ven* 
n=o n=o 
L 2r-2 . 
2 
+ exp (-2a°x) ) W aX 4 
n=o 


and the zero conditions on the boundary imply that for all r 
(1.50) Ula Ve ew 


Substitution of the expression (1.49) into (1.46) leads 
to the following equations for determination of the Nee in terms 


of the known constants Li; MeyoP 
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The corresponding equations for determination of the 


won are obtained from (1.51) by substitution of M41 for 


5 
Lig Orta for Berita » and 2a TOY ea a: 


The equations for determination of the ie are 
identical to the equations (1.34a) in which each nee is 


replaced by UL} the Us are determined from the relations 


(1.50). instead of (1.34b). In particular the value of Vio is 
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The next, and the improved, approximation to C(x,9) is 


C(x, 8) + MC, (x,0). 
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CHAPTER II 


SOLUTION IN CHEBYSHEV SERIES 


Ls Typical Characteristics of the Solution 


In this chapter, it is proposed to consider the solution 
of the simultaneous partial differential equations (1.1) and 


(1.4) in terms of. Chebyshev series. 


The reduced system (1.30) and (1.45) suggests that 
inherent instability exists in this case. The Complementary 


Function of the differential equation (1.30a), namely, 


g2c (ttl) 
eels) TS eel acitt)) = -gc{") Peach Us 
dx O fe) 
is of the form 
Ce ed) gue = Ae’ * Be ie ee 


However, the boundary conditions for x = ~ force it to the form 


etl) oa? 3 


(2«3) 6 e > 


va x 


eliminating the term Ae which, in fact, isa major contributor 













| Li #aTEAHD 


‘aaisue VaHaYaaHO MI MOTTULOR 





aotyules ait xebtenos ot besoqoxq st at ,zedqeds atds al... ©) 
bas (1.£) anoltaups Iatsnezeitibh Isittsg evoensti{umie sd3 30 
. -zoltzse vedaydedd to amiss at (@.) 


tad2 eteeggve (2),f) bos (O&.I) mosaye bsowbes sAT ale fo 
vsesasmsiqmod eit .ge65 ati ni edeixs ystitdasent dooxerint 


‘Yismen .(e0€.I) sotiaups Istaners?2ib ed3 to motsoaul 


nahh (is), S). tui? 
(% - f+). - o , 


mrot sd3 to et - 


si + 3 (S.S), 






mzo2 oif2 03 32 exx0l » = x ot eaotsibaoa yxabaved ads .sevewoH 
>. A | | il 





=. 


Los 


towards Cee’ 


for large values of x. Any approximation technique 
introduces this factor into the complementary solution in the 

form of round-off error. Thus the calculated solution may increase 
asymptotically in comparison with the true solution. In Chapter LV 
the problem has been considered in some detail and it is proved that 
if stability can be achieved, the solution is convergent. The problem 
of stability does not arise if there are no round-off and truncation 
errors. But this is very difficult to ensure for a computed solution. 
If a small error is introduced in moving from point A. to point 

B on the curve » (Fig 2.1) which represents true solution, so that 
the computed point is C, even if no error is introduced thereafter, 
the computed solution, for large interval, can lead to. N-- a 

result enormously different from the true solution. Even for small 
intervals, the final solution may be sufficiently removed from the 


true solution so.as to make it completely useless for all practical 


purposes. 





Fig. 2.1. Propagation of error. 





-€8 


suptaioed aoltemixoxqqe yaA «x ‘Qo esulev opzal x02 aia abyswod 
eg ci aoltuloe yrsinemsiqmos sd. ojat xo326% ainda assuborsal 
susetoal yam noltjvioe baisfvolss sdd eudT .to1rs tto-bavor to xo 
«VI tetqad2 al .motiuloe suri srt ditw noetzsqmos mf ulisstsosgayes 
Jedt bevorg af 32 bas [istsbh smoe al bezabhenos msed est meidorg afd 
mafdozq sil .3megtevnes eat notsuloe sd .bsveinos ed nso. ystlidawe tt 
notssonust bas 2o-~bnvox on sis exeds tt seizs tom ss0b ysttitdate to 
-tokiuloe betuqmos s xo} s1vens oF JluolIitb yisv al. etda aud .exzorre 
jnteq oF A intog moti gnivom ak bssubossat al zorxs [leme se 31 

tad3 os sotsuitos sutt ainsesiqe: doltdw (2.6 ali) | evtus of3 ao @ 
.t93%sexedd beowborinl et ae on ii asvo .D ak Jato. bojuq@os edz - 
6 -- M 03 bssl aso .lsvieint sgisi tot ,softuloe isasisgnlids eds 

Liswe rot mev% .nobsufoa eux sft mort gnststith ylevomrons sivees 
si3 movi bevomes yiinstotiive od vem motiuioe Isat? oft. .slavresat 


fsoisoszq [is toi easlseau yisasiqmoa 31 solem o4 es. os aoktuloe suxd 


. 29209 10g 


4 








24, 


The only way to meet such a contingency is to keep the 
increments in independent variables small. But in general this 
is also not a very satisfactory technique. It will be far more 
useful to maneuver the error into ripple-like character so that 
any divergence introduced at any step is negated immediately 
afterwards. Moreover, the fastest rate of convergence should be 
the target so-as to meet the delicate situation. It is with these 
factors in mind that the solution of differential equations (1.30) 


and (1.45) have been represented in terms of Chebyshev Series. 


Peas Chebyshev Polynomials 


The Chebyshev Polynomials belong to the general class of 


functions, the Ultraspherical Polynomials PO) (x) defined as 


nh 
(2.4a) PO (x) . c (ex?) se 


dx 


PES IESE Or 
where o is a constant and n the degree of the polynomial. 
They are orthogonal over [-1,1] with respect to the weight function 


eat CPAs oe 


For a= 0, they correspond to the Legendre Polynomials 


which are alternatively defined by the relation 





(2.4b) P (x) = 
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and satisfy the recurrence relation 


(2.5) 5 ai Oat perm el Daal erp lar eS: 
and the orthogonality relation 
uh 
(2.6) le P(x)P(x)dx = 55> (men), 
= 0 (m #n). 


An interesting property of the Legendre Polynomials 
concerning least square norm is stated in the theorem following 


the 


Definition 2.1. The Least Square Norm of a funetton f(x) 
with wetght funetiton w(x) over an interval [a,b] ts defined 


as 


b 
fio) || =f wle)f? (olde . 
a 
Theorem 2.1. Of all the monic polynomtals (t.e. wtth the 
coeffietent of oo” untty) defined over the interval [-1,1], the 
Legendre polynomial CP (ad where C ts chosen to make the 
coeffietent of a untty, have the least square norm, with the 


wetght funetton untty. 
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The proof of the above may be obtained from any standard 


text. (See, for example, Synder [23 ]). 


Another important property of all orthogonal polynomials, 


and very useful in the context of the problem, is stated as 


Theorem 2.2. If the wetght funetton w(x) does not change 
stgn tn the tnterval [a,b], the polynomial T(x) whteh ts a 
member of an orthogonal set of polynomials, possesses n distinct 


real zeroes all of whieh lte in the same interval [a,b] . 
For proof, see Synder [23]. 


The above theorem guarantees that P fx) has n real 
and distinct zeroes in [-1l,1]. Thus P ox) oscillates around 
zero with variable amplitude which increases as we move towards 
the end points. This corresponds to the case —- $< ee ae 
However, for. -1 <a < - 2, the amplitude of the oscillations 
decreases as we move from the origin towards the end points of 
the same interval. For a= - 5, the amplitude remains constant 
and the polynomials that correspond to this case are called 
Chebyshev Polynomials denoted by Tos) The corresponding weight 
function is 


l 
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which serves to keep the error small near the ends of the range. 


Se Sa | —  ? rn eee. |!) Eek, 
tot 7 7 





















198 


; 
i 1. 





busbnste yus most borisado sd vem svods sid: te°Rocrq oft 
| {EQ tebree potqmaxe wot ~962) | - sted 
cael 
,alabmoayloq Lanogodsx0 Ifp to vazeqorq Inet soqait ee 


as beistea at emoldoxg ed2 to gxeones edt “i ivteex = bas . 
‘re q 7, . 


spade ton esoh (nie nobtonnt tdprew sit Tl .S,8 mexoedt —~ _ 
pei foils (alin Ingmonylog sits ds] Sagreted oft ab apis 
sonttath «© sseasesog ,sinimonglog to iss Innogodine m Yo “tediom 


. [dao] Ja@catnd sme sft ab 935 Notts Yo Np seotes Inet 


.[ €§) tsbay2 sse .Rooxq xOF 


fsex a eed (x) ¢ sed? coatnatsuy mexood? svods SAT - 
bayors eodeiitoeo (x) I euiT .[i.1-] sb eeoxes dontietb bus 
ebiswod svom aw es ssasovoat dotdw ebusifqus sidatisy ditw oss 

- 7? > p> < - B85 9c3 - abnoqes 1109 etd .83ntoq bas edt 
asoisaiitose sit to sbujitlqms srs , =- >o>I- 62 < Tevewoll 

to etatog bne sit ebxawod migizo odt mos? svom sw e& asasetoeb = 
Jaademon entemex sbusifgqas odi , é - =-2 oF  avaedat emse od 
‘bells sxs sess ati 02 broqeerxo> ted? elstmonylog od3 bas 
adigkew gatbnoqeszyo2 sdT . (x) T yd betonsb elsimonylod veriey dod) 


60. nied o et mokzonut 


oe eS oo 









aes : 


a* “ w zo) i 


he. 


Thus Chebyshev Polynomials form an orthogonal set with 
respect to the weight function (2.8). Also they have all the 
zeros real, distinct, and lying in the interval [-1,1]. They 
oscillate around zero with constant amplitude unity throughout 
the interval [-1,1]. This is sometimes referred to as the 


"equal ripple property". 


A few standard results for Chebyshev Polynomials have 
been recorded without proof in Appendix A. The following 


definitions .and theorems are necessary here. 


Definition 2.2. The norm of the function f(x) tn the 


interval [a,b] ts defined as 


(2.9) |f(o) || = sup > Rec re 


(Qh ‘Ss Gy XS 


Theorem 2.3. Of all the monic polynomials of degree n 
on the interval [-1,1] , the polynomial PO eh has the 


smallest norm in the above sense and that 


(2.10) j27 7 (ad || = pia 


Synder [23] has proved it and stated it in a more useful form 


as: 


Among all the polynomials of degree n, wtth maximum 
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norm untty in the interval [-1,1], T fa) has the largest 


leading coeffietent, namely, peed 


This immediately leads to the 
Theorem 2.4. In the class of Ultraspherical Polynomials, 


Chebyshev Polynomials yteld expanstons which display the strongest 


posstble convergence. 


For proof, see Synder [23]. He has further proved that 
the Taylor Series which is a limiting case of Ultraspherical 


Polynomails for a=° display the worst possible convergence. 


2p LruncaLLonshrror 


If a function f(x) is expanded in terms of Ultraspherical 


Polynomials, that is 


(bah) f(x) = } GoBack) 
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and is approximated by 
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the truncation error is 
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(a) 
C2EA13)) € 6s) = CP. (x) 


For Taylor Series, a=© and € 6X) = Cxy which 
increases continuously as we move from the origin towards the 
end points. For Legendre Polynomials, a= 0 and € 6%) = CP fe) 
which is lowest at the origin but increases as we move towards the 
end points + 1 and-is of maximum magnitude Ci For Chebyshev 


Series, a= - 


Nile 


and € (x) a CTX) which has an equal ripple 
cosine curve. . (Fig. 2.2). Thus error is more evenly distributed 


over the entire range. 


(Error) 





Tilo ewe. 2 meeedqual ripples curve: 


A Taylor Series expansion of a function will have large 
coefficients even though the size of the function throughout the 
range may be small whereas the Chebyshev form will have much 


smaller coefficients and hence a greater facility for significant 
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. 1 = Ge) ys 
| | eal 
dotdw "x9 2 (x) a bos * = © vestise sofysT xoF | 

ody abyawod atgtro edz mori svom sw es ylevountsn0o equneange 

(x) 4,9 = (x) “? bas O=. ,alatmonylod exbasgel roi .eintog bas 
ed3 sbrawod svom ow es ssasetoni tud atgtyo siz 35 deswol et dotdw 
vordayded) sof =. 9 sbustagsm mumixem to ef bos I + asntoy bas 
efqqtt Ilsupe as sad doidw (x) 79 = (x)_3 bas = -= Pa .asize? 
besudixaetb eaneve stom ef torts audT .(S&.§ ght) -ovwus smtaoo 
.930Bt stiios sda seve 


SO). 


truncation with little error. 


Any finite range [a,b] may be linearly transformed 
to [-1,1] to which pertain the Chebyshev Polynomials. Similarly 
to the range [0,1] pertain the (shifted) Chebyshev Polynomials, 


defined as 
ie TY a Ge) A Oreess) = a CR) A ee 

; n n nes. ea a ae =* i> 
with corresponding weight function 


il 1 
ae: — 5 


(2.15) egbe [12x -)) = fenton) 

It is interesting to note that for the approximation of a general 

function f(x), we should expect, with the same number of terms, 
* 

to get a smaller minimax error with T (x) tne [10.1] 0) thanewith 


TG) in the larger range [-1,1]. 


2.4. Fourier Series 


Since Fourier Series are based on use of orthogonal 
functions and are intimately connected with Chebyshev Series, in 
this section, it is proposed to consider the same and compare the 


convergence properties of the two. 


The trigonometric functions 1,cos x,cos 2X vee ey oun aks 


sin 2x,... are orthogonal with respect to the weight function 
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unity in the interval [-1,7].. Hence the trigonometric 
polynomial 
n 


‘lh 
(216) P bX) = as + Es (a, cos kx +b, sin kx ee 


with the Fourier Coefficients defined by 


T 7 
Cran) s =f f(x) cos kxdx ,  b ==) f(x) sin kxdx , 


—T Th 


is a least square approximation to f(x) with unit weight function, 


in the interval. [-7,7]. 


At any interior point x of the interval in which £(x) 
is bounded and has a finite number of Maxima and Minima and a 
finite number of non-coincident discontinuities, the Fourier 
Series converges to = [£(x,)+£_)1 which reduces to f(x) at 
a point of continuity. However, the rate of convergence of the 
Fourier Series, which in other words is the rate of decrease 
of its coefficients, depends on the degree of smoothness of the 
function, measured by the order of the derivative which first 
becomes discontinuous at any point in the closed interval 
[=t,mhe EWen if f(x) is perfectly smooth, it may have 
terminal discontinuities with regard to the Fourier theory which 


in turn affects the rate of convergence of Fourier Series. 


Consider the Fourier expansion of the function rix) 
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defined in the interval [-1,1] and expressed as the sum of even 
and odd functions £, (x) = = [£(x)+£(-x) and £, (x) =5 [£(x)-£(-x) ] 


respectively. Thus 


1 4 a oy! 
£ (x) 5 a + ) a, cos kx, a, = = f(x) cos kxdx , 
. kel fe) 
(2.18) 
At) we Ce = Sa rl f(x) sin kxdx . 


All terms of the sine series vanish at the middle point 
x = 0 where £, (x) = 0 and at the terminal point x =T where 
£. (x) = = [£(n)-£(-n)]. Hence unless f(17) = f£(-7), the series 
will not converge at the end points of the range and will 
converge very slowly at the intermediate points. Of course, the 
first derivative of the cosine series vanishes at the terminal 
points. But this is far less serious than the discontinuity 
in the function itself and we should expect the cosine series to 
converge much faster than the sine series. To be more specific, 


on integrating by parts, we have 


TT TT 
(2.19a) f f(x) cos kxdx 2 [f£'(x) cos kx]? ~ +] £"(x).ecos kxdx , 
fo) k k fo) 


T TT 

1 ik 

(2.19b) f £(x) sin kxdx = - 7 [f(x) cos kx] ¢ as J £'(x) cos kxdx . 
(e) 
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For large values of k the integrals are likely to be small, 
being dominated by the oscillations in cos kx. Thus ultimately 


the cosine series converges like Sat and sine series like we 


Considering the range [-1,1] and making the transformation 


xX = cos -6, so that 
(2520) E(x) = f€toes 9) ie 2(G) mM MOS 9 < 7 


The function g(6@) is even and "genuinely" periodic with a period 
2m. In addition, if f(x) has a large number of finite derivatives 
i (-1,1), so has ¢(@) in [0,1]. iInterpretting the Fourier 


cosine series 


(2.20) g(6) -5 a + ) a, cos ké@ , a. 


T 
= i g(8@) cos kédée , 
fo) 


in terms of x, we produce the Chebyshev Series 


1 
2,-% 


(2.22) f(x) = tat J at), ast f (ex) *£G)T, (xddx . 


This has the same convergence properties as the Fourier Series 
with the additional advantage that the terminal discontinuities 
are eliminated. Also for sufficiently smooth functions, ay has 
the order of magnitude oem ace) to which eee stands no 


comparison for sufficiently large values of k. 
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2.5. Semi-Discrete Solution for C (x58) and T (x58) 


Before proceeding with the solution of the system, we 


wish to prove the 


Lemma 2.1. Jf 


-ke ,1 * 
fed ar "ta ,P ite. ta 2) dee 
(2.23) 
= ot Pen perth Vie 
Guo gua nn 
we have 
(2.24) 





Proof: Differentiating both sides of (2.23) with respect 
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where 


' _ ' = = ' = ’ = 
py bel 4rb zee LCL) J and be b nF. 


ntl 
Comparing the coefficients, 
a= Kb. + bi) > xem O(1l)n& 
Therefore, 
rales ortlon ASS og gy HS s 
i.e 
_ 4(r+1) Ae Bt? 
Meo Ee a a AP ey) kk? 
Yon nd, oe 60S oath = b+? = 0; atl = a+? 0 


This proves the Lemma. 


Reverting to the solution of (1.30), the substitution 
x = Ry. where R is the depth at which concentration and 


temperature are required, the same are modified to 


qz2c (eth) 


fh 
(2.25) +25 -aci* ie SPC Mie 6 
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ger (etl) 
(2.26) + 2 yrer) ps yn?) " cet 
R dy 


(x) (r) 
where Co = Co (y58,)> TS = T (v8): 


The boundary conditions (1.32) are. modified to 


(x) (r) = 
(2.27a) Och x AO Sree 
(2.27b) c6 (y) =i05, Tay) =O, 
(2.27¢) ieee oA TR EER Saye 
Roo ° R00 ° 


For r=0 and 1, the equation (2.25) yields 


(2.28a) eng) = £,:exp(-R va y) > 
(2.28b) c8) Gy) = rales |) GD ae 
fe) oVa 2 


For general values of r, the solution is of the form 


(x) cet * 
(2229) Ce (y) = exp(-R Va y) 2 a a 5 


where prime in the sigma implies that first term is to be taken 


by a factor of half. 


Multiplying both sides of (2.25) by sexpt—k va y), the 
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differential equation is made exact and we have 


ac (rt) 


(ee Sa Cai exp(-R va y) = -8 f eva Bey 


(2530) © [ ly 


Grad 
= = ' * 
exp(-2R Yay) i Beery) 4 


where the use of Lemma 2.1 and equation (2.29) gives 


c 2(n+1) iy 


oS ren eee) r,ntl i r,nt2 D aR 


emit. t 2, wees A. =O) PLOL es rie. 


The constant of integration vanishes by taking the limits 


asueRae i =. 


Equation (2.30) may further be integrated to give 
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C22) exp (R Vo y)Co 
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=R f L AL Toy 3 
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peagemeys = 0 in (2.32). Also recall that Gia (0 = f. 


+1 
whence 
3 n-l 
S10 77h ) erin 2 ]. 
n=1 
Combining (2.31) and (2.33) gives 
mee Raed Pe a earl! 
ay eae 
; vo, ; 8n vo 
Changing n to n+ 2, subtracting and utilising (2.33), 
we get 
E _ 2(ntl) hs x agpeliaasesl 
eared hte Rays eal h reall rtl jnt2 
(25.34) 
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Now to determine T Cy 98) for xr = 0, the equation 


(2.26) yields 
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For general values of r, the solution OL ec eeG) oeleaot 
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r-l 


C 
7 es 3 * : : ¥ 
y exp(-2R vy y) ) Baaty + exp(-Rvy) ) Seiya A 
n=o n=o 
where 
23 = — 2 = 
( 2) eerie Yer on aa 2 O(1)r 
Integrating (2.38), 
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where taking limits, as R-+~, the constant of integration again 


vanishes. Also the use of Lemma 2.1 yields 


mee (ira) r,n 4r,nt2 
fay, oven 1, R Petal es) ala 2R # 
d -d 
peed (ort) eect near 
(2,415) oie “QieveR rtl1,nt1 ss veeiear2 vR 1 
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that is, 
(r+1) . 
RY = 1 
C2242) To = exp(-R vy y) L ete _ 
Te? § i 
= ' 
+ exp(-R va y) ) Seine ne 
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In. (2-42) set y= 0. “Also recall that To (y=0) = Bt 
whence 
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The constant of integration vanishes because each term in 
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Again multiplying each side of (2.54) by exp(2R Va y) 
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CHAPTER III 


SIMPLIFICATION OF SOLUTION IN CHEBYSHEV SERIES 


3.1.9° General 


In the last chapter, recursion relations have been 
derived to give the coefficients associated with the next-step 
differential equation in terms of the same of the previous-step 
solution. In the present chapter, these relations have been 
expressed in a much simpler and direct manner by employing 
matrices. The next-step solution depends upon the derivation 
of the inverse of next-order matrix from the previous-step matrix. 
The matrices are lower triangular and therefore it is not 
difficult to find their inverses by the Gauss Elimination method. 
However, they are generated in a particular manner which makes it 
convenient to derive the inverse of (r+l)fth order matrix from 
that of rth order matrix. Thus stepping to the next stage 
amounts to extension of these matrices by the addition of last 
row, the last column except for the diagonal element in each case 
being that of zeroes. Although the present technique is much 
simpler for automatic programming on a digital computer, it will 
be observed that it is not as fast for a desk calculator for which 


the relations of the last chapter are more direct and easier to 


handle. 
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Any apparent difference between the results of the last 
chapter and this chapter is only superficial. In the last chapter, 
the (r+l)th step coefficients were expressed in backward recursion 
relation in terms of the (r+l)th step coefficients which follow 
and the rth step coefficients which were known. If the recursion 
relations are conclusively employed so that (rt+1l)th step 
coefficients are expressed in terms of the rth step coefficients 
only,we have the relations of the present chapter. Thus, in case 
there is no round-off error, the results of the last and present 


chapter must agree. 


3.2. Semi-Discrete Solution for C (x58) and T (x58) 
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and a similar result may be found ciabenGy) by replacing r by 
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Substituting in (3.24), and comparing the coefficients 
x 
of exponents followed by comparison of coefficients of like 4 DAY 


we get the relations 











2p vee’, =i i i 
(3.26a) ~ R oe eS + 2 Cee Bul 5 (Cn O(1)r-1) , 
Ose © (een ae — eee AR Ly 
; u Y rtl,n R rtl,n Re rtl,n ron relia 
(n’ 2-0 G)2r)*4 
4 va, os 
a) SEM og ke ae tae Wrtl,n R- “rtl,n ae % Mr+1yn 4 
kn ="0(1)2r)43 
where 
(3zbd) na =O) (forjen > xr; Le = 0 = Fe foremarecra— we 
Again using (A-23) and (A-24), we get 
PS 
(3.27a) R AL Got ee) orbit OP a ete pa ote es 











wd 2 eatin % cy 5 bauot: scene zettate » bas 
evods sd¥ mt ca 


5 

‘etnetoliteos 4d? gatxsqmoo bus ,(#S8.€) at amtausiiedve y 
.e' T sit Yo eanstotiieos to doalzaqmos yd bewollot esmenoqée 20 | iy 
enotialex eda jeg Sw 


tp oe aS | | 
« (ia)0 = a) 88" = tay Sy + kee” 5 (sa .£) . 
‘age * aa’ " an’ 5,7 tn” & rag OFOW SHY) anne) 
; CxS(L)0 = mw) : a 
= thes i py 2 4 
a a a os = vee ca ater (2d8.€) 
any 


« (s8(1)0 = a) 


oe xo? MO: ie ag I Ry OOM ew u 


ft, 
; rt .* S y ; @ a P : 
ye we: 4 
aN ; ra pais ; 


is a” 






a 


ieee 


pial 





OZe 


+28 cary Cntv2) u +2 (n+2) (nt4)u 
R 


rt1,nt+2 rt+1,nt+4 
+ PDN aay, yt = Tees 
(nem Cl) re -al)ae. 
(3.2/0) (y+2 YaY)V er on - = CS ..§ CEG 


a8 Cae ie) Be ay 
R 


r+1,n+2 rt+1,n+4 


+ 3(nt3) (n+6)v ] = Reve 


eee + 
erie ator n rien : 


(ni nO Ci 23) 


4 va 
(3.27¢) 30W 41 on - rnd Ga On CC CD 
- a" (hel) (ata) +2(nt2) (nt4)w 
R- r+1,n+2 rt+1,n+4 


fe 3(nt3) (nt6)Wiiy neet est! = hal ap Ment ve © 
(n = 0(1)2r) 


Paton sO el, 25.0) eto. seu 












ican . a, ; on oi tk ‘+ 
sata degC) (St) St, pu (Se) CL) % : 
tet Lee typ ite Orn Eee + 
ato Te 
: a ~ 2(1)0 = a) ‘ee 7 
[ante prep MCCtOt ppg Cltn) 10 wee mn) 
7 | al 
debi, fa PM) CSOD 14g MES IG) htm) LI] ae . 
© ggea® * aya & b+ +o Magee: tO) COPE + - 


« (s$(D0 = g) 43.00 a 


| ov d 
Coe Hare cegelCtite sp wernt py ywilem) 1) =e * a reg¥OE CONSE) 





lias pea? OO) (S+a) St. 


ba, tea? S42) (L+n)£) = + 






- ~ «¢@ atts” ta," Les ety gt py COtn) (+a) + 


_ 
-— 


- 


. 





omy 39GS@)0 = a) LE Wel yahen : 
1’ ee —_ ; , 






63. 





2 
3.28 Mee oe 
( a) C2 Shas rs 16 Ns ey 0 tha aay) 4 
g Re 
3.28b eee = ee eee 
( ) rboleried “Wav re le oreeor1O ee \Yrso’e. tee Me, 27-20? 0) 
- Re (2 Q Q ) 
16% rFlso rtTeler!” @r+1 27° ° 
BR? 
(3.28) Seo een al aaa Uerior rele Sie OR yaphy 00 bs, Gee 
_R (m m m ) 
16 ‘ rtl,o rt1,1°°° rt+l1,2r 
where 
G5 0 O 0 0 0 ee 
E19 B 5 0 0 0 Om. 
-1.1.2 2€45 b 5 O 0 O . 
(3.29) Ertl = 3€ 19 =1L. 203 3E 15 G5 0 OF .3. 
-2.2.4 449 -1.3.4 4E 45 G5 0 ee 
Sage LiL Ge BE Ry eee dan ec 
eee 12 ae 12 DE Nats 
2r+1x2r+1 
with 
pw wit? Yay) R are 
(3.30) D 16 Pare ay 2) 2 





ed 


aie = 
e Creag een Wy gh = ake Sy rae” aae ee eer 


s 


a 
He ccgt.a¥***£<3%0,2" oe 


r Cog, tta® weer 


(0. Ores. -. the * ‘“ 


. € 








3 fe "tayo 


an a 


&.€.I- aA Os A & SS —] . 


see <3 ytd ¢.A.L- 


Si 


das aces” °° Lytteo ite” ~ 
g %) Sg ~ 
t beso. ite™ OI 
ws ea 8" peo tee 
$ 
A 
i) BIL - 


3€ €.°,f[- 3€ 


St 
Si 


ae, C,658-. SE 















64. 


Se 0 0 OE aie) 4 ae 
2e4 53 0 0 0 a 
=-1 0,2 4e 4 b3 0 0 ots 
(Gerry : 
2r+1 be, -1.2.3 be) b3 0 6 
2rtix2rel 
with 
ED 
(3932) é, 7 =30R 
16 
The technique for deriving the matrices ED etl and Foe 
is similar to that for deriving Coat! Forel may be obtained 


directly from Eo e+] by replacing E5 by E3 and Eo by Ey: 


Thus we have 


= BR -1 

(3-338) g{U iy abpbeoeyed Visi 2ld eine mblesots, ave: Us feitie > 

af 
(3.33b) LV ey onted 1 * BS Yael qeobeherolrta a Mepztazne: © 

- (2 Q cae Vine 
rtl,o rt1,1 r+l1,2r 2rt+l1 ’” 

Re 

(3:33e) gy Oia dt eM eet, 2e) 7 16 BO, Meas Mr, 2r-2°920? 


















a) 
1 0 0 0 °3 1s 
wae QO. 0 ee ios $i£.4- " ne a 
‘ 3 q . 
: *** 0 e? ~™ £.8.i- + f S j 
eee e° 1 @,€.I- 1% ».£.S- 
featutans *.* '* . . . . * . . * . 7 
djiw 
“ave- 
ae. 3 (S€.€) 
al ce ae 
t aa? base L 4ee sootijam of3 gnuiviseb 10% suptodoss sAT 


beakesdo ed yom 5.01 «)..9 gatvizeb soi tadz o3 selimte at 7 
+3 yd ¢? baa ¢? xd <3 gatoalqss yd tere mozit ylsosakb 


eved ow eusT ei a / 
| | a 
P| ao oe - 
fa 3” eae en Ped =e On eon teat tee? (at&.€) 

sets ~ 


7 
7 
8 


s a 
tt aad ‘aa 4) 5t * Coc to” ooe .tex"0,.d42” (dEE.E) j 


7 - *; 





—- 


a 














sa ok 
eats 


ieele sean wa 








OD. 


a 
oir peor ore sarki 


= 
Recall that Creek (0) = 0. Hence from (3.25) with 


yr replaced by r+ 1, we have 


C3 950) ""u = 2[(u ore ( 


f a) IPA GC 
Er. .o rer “rtl,2 aang ee 2 Thentiy poh 


+ ( 4. Jeti (w 


Mera) Win Sh) rel, rt1,2¢° 92)! 


The equations (3.33) determine gy (5) completely. 


Pinally. as? in (2259) 5 


ir p-1L 2p-2 
= = a i] nel ! 
(3534) C1 Cy by) exp(-R va) ) aan + exp(-Rv) ) ae 
n=o n=o 
2p-2 
+ exp(-2R Va) - Ww : 
sey ee 


3.4. A Note on Inversion of Triangular Matrices 


It may be observed that all the matrices involved in the 
previous sections are lower traingular. Moreover in each class, 
a matrix of any order is systematically connected to a matrix of 
next higher order. Hence the following technique for inversion of 


such matrices should be useful. 


Let the inverse of 















= 


= | 7 
i- i ‘ 
* pase peg? geno tee 


datw. €S.8) g@eSPennsH . 0 = (0) ati dais ttaces 
sved ow i +s yd beoalger + 


I ome 2 v. 
Couaaogtea? E> Grate ts tase sta T® * ota” OORE 
- 
+ The race * Cte tee ee” * 


.visssiqmoa (yx) a ‘gaimysjeb (f£€.£) amoksaupe ediT 
(22.8) al op vettantt Eo . 


v "7 (vi-)qxe + Fe ge BY f-)qxe = (Ley) 9 (dE.€) : 


( 
Org 
S=qS 
' ad? |S Co\ @S-)qxe + . i 


to molataval mo stow A LAE 


ail at beviovst ssptztam srs Ihe 2ed3 bevieedo 4d yeu 3I | sz 
+sealo doses mk reveox0M .1elugeias3 rewol sxs enolase evotvexq : 
20 xbzgem 8 02 bstosmnos vilsotiamstaye et z9b30 yaa to — 
Re eetamat tot suptmions. gntvoliod 222 socer .tabro xsdgid usa 
Sunes, yas tina we:Aiciats. seni 






. 









Fo 0 
a a e e 
(3.35) Korps 
= 
rl r2 - ay or 
be called a - We can partition 
| 
A. Oe 
= =S& 6slac aae 
(3.36) Ans ! ’ 
Le eas Art rt 
where 
(3.37) Bo nee eel 8? 


Let 

- om 

-l 

(738) A 4d ; 

ty resale 
so that 

1 es Cee 
(3.39) Coe Gere) 
LB te ely Oher+pekh etd 


66. 


it 






3a 
oa 
Poa . 
gotzixeq map sW .7-A belles od 
‘ 5 . 
tx” =” } 
-_——-— - — = Me —_— = ° 
i+a* 
i 
rey +47 
re gytex*t ce” ° 47 
joel 
a. e,.3 
(8€. 
7 i (Rg s 
e penier 





> yas Ye cttign e 


> ae 


D eee att a ee (wales sangee 


iin ue rca 


7 ae | 






a 







Lt 
At 


67. 


Hence 
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3.5. Economisation 


The conclusions derived so far regarding Co: qT) and 
C, are valid over very short intervals of time. The technique 


uses the known values along Ox and at A, (Fig. 3.1) to yield 


1 
their values for all points on A,B, . In the next step, A,B, 
is taken as a new x-axis which, therefore, with the value at A, 


yields their values for all points on A Proceeding 


985: 
similarly, finally we get the values for all points on AB, : 

Thus, provided the error does not accumulate to produce divergence, 
the process may be applied over a relatively larger time interval. 
But as suggested in section 2.1,Chapter II, the time increment Aé 
must be kept small to avoid instability. This, however, 

increases the number n of time steps. To arrive at the final 
solution, the matrices A,B,C,E,F of order n,n,ntl,2n+l1,2ntl 
respectively are generated. This in many cases may prove to be 

time consuming, if not impossible, keeping in mind the limited 
memory space even in the largest computers. Moreover, it being 

a property of the Chebyshev Polynomials -- and our preference of 
Chebyshev polynomials over Taylor polynomials was motivated by 

this property -- that in a convergent series expansion of a function, 
the latter coefficients tend to decrease at the fastest possible 
rate. Thus it may happen that after a few steps, the later 
coefficients are less than our choice of the criterion and hence 
negligible for all practical purposes. It suggests that some 


"economisation'"’ in the process at an intermediate stage should be 
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a desirable feature. 


Let us suppose that for some r, the calculations yield 
a3.465)™ "zero" (i.e. less than the criterion). In obtaining the 
3 
next set of coefficients, we get 
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Stir = lore, ee yr-l ' 





(3.43) 
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Provided < 1 (and once so, it will be true for 


all subsequent steps), we have 


(3.44) EN ate Oe Zero. 
This sets up the chain process 
(3.45) "zero" = 


Art2 r+ z Frt3 rt2 a 


Thus the (r+l)th step onwards significant coefficients are given 
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Fig. 3.1. General Scheme for the Solution of the System. 
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and so on. Thus the reduction of last coefficient a. py at 
rth step, to ''zero' implies that the dimension of A. be not 


increased before switching to Grplpithewscep derivation. 
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at some rth time step, the dimension of none of the matrices 


A,B,C,E and F may be increased in shifting from the rth 
step to (r+l)th step. In other words, after completion of 
each step, it will be relevent to examine if the above conditions 
are satisfied. If so, the dimension of none of the matrices need 


be raised. If not, proceed as usual. 
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CHAPTER IV 


CONVERGENCE AND STABILITY 


4.1. General Terms 


In this Chapter the Convergence and Stability properties 
of the partial differential equations (1.25a), (1.25b) and 
(1.26a) will be considered. Hence it will be appropriate to 


have a general view of the terms involved. 


In the solution of a partial differential equation in 
which the dependent variable U is a function of the independent 
variables x and t, the partial differential equation is 
approximated by a difference equation or an ordinary differential 
equation. Let u denote the solution of the approximating 
difference or differential equation. If |U-u] , often referred 
to as the "discretisation error'', tends to zero as the increments 
in independent variables tend to zero, the solution is said to 


be "Convergent". 


However, the round-off errors may affect the calculated 
solution u to produce the numerical solution N of the 


approximating equation. This gives rise to two possibilities: 


(i) The error introduced at any step may increase exponentially 
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over the subsequent steps. Thus even if no error is introduced 
at a latter step, the procedure has the tendency to magnify the 
error at each step. Consequently, the numerical solution N 

may differ exponentially from the solution u. If by proper 
adjustment of increments in independent variables (which places 
limits on the values they may assume for producing a useful 
solution) it is possible to avoid this type of adverse situation, 
the method is said to be "Partially Stable". But if this is not 


possible, the method is called "Inherently Unstable". 


(ii) If the error introduced at any step tends to decay in 


subsequent steps, the method is said to be "Stable". 


Thus whereas the convergence implies the decay of error 
| U-u| , the stability assures the decay of the rounding error 
|u-N|... In other words, though each affects the true solution 
U, they are independent of each other. A process, in which 
calculations are performed uptominfinite number of decimal 
places, or in which no round-off errors are introduced, is 


necessarily stable. 


Convergence, in general, is more difficult to investigate 
as it requires knowledge of the limits on the derivatives of U 
which, in general, are unknown. However, it is more relevant to 
the solution of the differential equation. Co-existence with 
instability is possible under certain circumstances, but 


convergence is a necessary condition for the usefulness of any 
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scheme. In some cases the round-off errors may persist linearly 

thus destroying a few final decimal digits, yet giving a correct 
solution upto a few initial decimal places. In case a scheme is 
unstable but convergent, it is called "Weakly Unstable" as 

compared to the case of "Strong Unstability" which implies non- 
convergence as well as instability. The former is unstable in our 
sense for some types of differential systems whereas the latter is 
always catastrophic. In case of weak instability, if we are not 
covering too large a range or if we are keeping a few guarding digits, 
a useful solution can be produced. Stronginstability does not hold 


any such promise. 


In the next sections we shall discuss the convergence and 


stability of the solutions Co: TS and Cc, of partial differential 


equations (1.25a), (1.25b) and (1.26a) respectively. 


4,2. Convergence of C (x58) 


The approximation (1.29a), viz. 
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transforms the equation (1.25a) to the form (1.30a). Thus the 


true solution of the approximating equation 
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at a depth x and at time (r+l1)A6 is ott) = c (x, rt+l* de) 
whereas the true solution of the partial differential equation 
Gi 258). viz. 


9°C_ (x ,xt1+A8) 


1 C(x, rtl-Ae) 


(453) Sag Ta) aE Of aa Ca LC (x, rtl*A6@) = 0, 

ox 
is ones = C(x rtl- Ae) . L£ we denote the error at the 
depth 9x fat time r+ 1°°A0 by. E47 (*) » or simply by Ee? 
we have 

Le rrl) ae rs 1) 
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For a particular r, r*AO may be regarded constant and hence 
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2 
die aC Ste 
(457) etl oF + BE _, He ee Ta 
dx” rt ny D ele) 


C(x, rtl* Ae) - C(x, 4-A8) 
A® 





ao tr et an 
D 08 D 


aC (x, rte »A8) 
06 


Z 
Ae ) C (xs rte -A@) 


= —=______ 
D 992 


where §0.< ¢ < «!S2 1. 


E. vanishes at the boundaries x = 0 and x == for all non- 
negative integral r. Further, since the term 

9°C (x rte +8) /367 denotes the rate of acceleration as far as 
the increase in oe in the @-direction at a specified depth x 
is concerned, it is bound to be quite small and almost constant. 
Therefore,let 
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The error giving equation in the final form takes the shape 
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For small values of x, the first part is approximately 


2 -Va x 
(BL tXB-X Be 
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Hence, provided the increment in x direction is small 


and of the order (4.15), the solution converges. 


Weoe) OCabLlity of C (x58) 


The differential equation (1.30a) may be expressed as 
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(4.21) (Da) *e(*) era: 


Sait 
the solution for which is 
(4.22) @(t) = (A +A,xt...t¢A _x Je 
fe) o 61 
+ (BOt+By xt. . ee 


which shows that error tends to increase exponentially for any 


x. In other words, the solution is unstable. 


4.4. Convergence of TS (%,8) 


Here the partial differential equation (1.25b) , viz. 
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Defi 
efining the error 2 coal as 


ect) (r+1) 
(4,26) Ba esd 2 iy. : 


substituting in (4.25), and simplifying as in Section 4.2, 
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Initial and boundary conditions on E. are as for Cy with a 


and 8 both replaced by y . Hence,in general, 


o 2 r -Vo x 
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=) = Gam) a 
where 
(4.28b) X =x vy/2, 
4.28 B, = eS) 
(4. 28c) ie 


Hence the first part vanishes for x ® 3.236/vy whereas the second 


part approaches zero in the limit as A@ + 0. This establishes 
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the convergence of Ty ; 


Heo. otability of T (x58) 
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Let be the calculated value of ra 


Pee as obtained 


on solving the differential equation (4.25). Therefore, 
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Subtracting it from (4.25), we have the error giving equation 


(4.30) — - vee.) #0, 


which is the same as in the case of Cy of Section 4.3 with 


a and 8 replaced by y. Hence the solution is unstable. 


4.6. Convergence of C, (x, 8) 
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Hence the solution is convergent. 


nie baDLLity 2Ot C, (x56) 


Here the true solution of the approximating equation 


satisfies the differential equation 
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which is same as in case of Co: Hence the solution is unstable. 
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CHAPTER V 


CONCLUSIONS AND SUGGESTIONS FOR FUTURE RESEARCH 


5.L. ~Conelusions 


To test the validity of the schemes expounded in 
Chapters I (where the solution has been expressed as a Taylor 
Series), II and III (where the solution has been expressed as 
a (renvanave Series), programs GAS1, GAS2 and GAS3, recorded 
in Appendix B, have been written and tested on the IBM System 
360/67 at the University of Alberta, Edmonton. The numerical 
results have been computed for absorption of Carbon Dioxide 
in 1M NaOH at a pressure of one atmosphere. It is believed 
that the realistic values of the constants involved are as 


follows: 


a = 7000 
-3 

a, = 7 29 2eece LO 
a, = 2.9 x re 
Ty = 520. F (Initial Temperature) 
Ep eke tateci ieplbienole/fte 

a 2 
es a.r.A0 F 
2 = 4000 
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8 fies iceee (Mass Diffusivity) 


Dewi 2e25exnlOe 
A = 71280 Btu/lb. mole 

b = 9.45 x ome Btu/ft. sec. °F (Thermal Conductivity) 
p = 62.39999905 1b./ft.° (Molal Density) 


C c= L460) weBtv Alb aeer (Molal Specifie Hest) 


If time-increment A@ = ere it follows that 


oi ee0w44622200—e 1Oht 


B = 0.4444444 x 1077, 


ROG 03nd 0 aa 


A = 0.30171432 x Mee 


For the convergence of C, it is required that the distance-step 


size 


But to obtain a convergent solution for To» the distance-step 


size 


Lee Eye Eas < aig 


which is ten times as for C. Thus, the distance-step size has 


to be adjusted according as the convergence of C or T is 
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desired. Their relative magnitudes depend upon the system under 
study. Thus in many cases, as in the present one, the simultaneous 
convergence of C and T may not be possible due to different 
distance step-sizes required for each. In the present case, the 
convergence of C was desired and it was observed that convergence 
with eretyeney” Series is much faster than with Taylor Series. 

In the case of Taylor Series, it was required to make 3.55x10 ~ fo « 
correction in On to get a better approximation ve + Mx C.> 
whereas similar correction for the cHebyatiey Series ae DSexace te 
Thus in cases where it may be incumbent to calculate Cc, or even 

C., with Taylor Series approach, the value of Cy calculated with 


3 


* 
Chebyshev Series approach may give equally reliable results. 


The recursion relations of Chapter II involve many 
multiplications at each step. This has the tendency to increase the 
rounding error, eroding with it the relative advantages of 
ano Series over Taylor Series. Hence, in practice, the 
method expounded in Chapter III proves to be easiest and safest. 
Whereas Sullivan's [23] method is applicable for only one milli- 


second, this method produces a convergent solution upto ten milli- 


seconds. 


Is she a Suggestions for Future Research 


The results of Chapter IV lead us to the conclusion that 


system (1.30), (1.45) is convergent but unstable. This requires 
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a great deal of delicacy in approach in dealing with the various 
elements of the solution. It would be desirable to evolve 

a strongly stable scheme so as to achieve fast and simultaneous 
convergence of C and T with consistent limits on the values 

of Ax and A®@. Whereas the speed of convergence may be 
accelerated by ©MD Loying Chehyaheu) Series instead of Taylor Series 
(if at all such a usage becomes necessary), the other objective 

can be achieved by trying entirely a different type of solution. 
Unfortunately, Lanczos t-method and the method of finite differences 
are not suited to the solution of system (1.1) and (1.4). Lanzos 
t-method gives the true solution as a finite series of a slightly 
perturbed ordinary differential equation when an approximate 
solution in finite series of the true differential equation 

is not possible. The system (1.1), (1.4) does not qualify for 
this method and true finite-series solutions of (1.30), (1.45) 

do exist. The finite-difference approach to (1.1) and (1.4) 


is limited by the fact that the liquid extends to infinity. 


The most promising approach might, therefore, be the 
modification of the semi-discrete method in which an alternative 
to approximation-equation (1.29) is used so as to produce a set of 
ordinary differential equations which are stable and converge 


rapidly. 
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23. Synder, Martin Avery. Chebyshev Methods in Numerical 
Approximation. Prentice-Hall, Inc. Englewood Cliffs, 
N.J. - 1966. 





APPENDIX A 


STANDARD RESULTS FOR CHEBYSHEV POLYNOMIALS 


In this Appendix, a few standard results for the 
x 
Chebyshev Polynomials dealing with the range [0,1] -- 
in Chapter I-IV and in general referred to as Chebyshev 


Polynomials -- have been recorded without proof. 


The general Chebyshev Function is defined by the 


relation 
ey T, (2) =5 (owe ony aa ey ae 


For i} =n , the nth Chebyshev Polynomial associated 


with this symbol is 
(A-2) T(z) => [Cet Y22-1)" + (2- Ve2=1)"] 


Substitution z*#= x (real) and the change x = cos @ reduces it 


to the form 
(A-3) T OX) mpCOB@NO XMM eCOg 0.5 .(—lecex cel) 


For the range [0,1], the change of x to 2x-l gives the 
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bststooses [etmonylod vedeydedd din sdt , a= k 20% 


ak lodmye eid3 d3iw 
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* 
Chebyshev Polynomial 


& 
(A-4) T(x) = T\(2x-1) = 21° (Vx) -1* 7, (/x) 
= cos (nd) , ¢ = ee ea yf) Ss cs aay 
Hence 
(A-5) DG)T G@) = 5 (IT, Get GO), 


which for m= 1, gives the standard recursion formula 
6 i 2(2x-1)T (x) + T. 0 
(A-6) Tyee) = 2(2xr LT Gt T(z) sO 


From trigonomety 


Ti Gena ; 
* 
(A-7) T, (x) = 2x-l ; 
T” (x) ee CP) P? (1x) “% 2B) x? (1x)? aw 


The integration formulae 


GN) gt = 1 - Gin = (x8) GT 


2 %20) . UHxS)*7200 = 9 «. (om) 200 = 


© (G0), HO), “T) zs (x) "a¢x) ‘ 


aluoso2? moteiwos1 bisbasze sft esvig .t = a x03 


1% 
Om Gd, OT + (x) T(L-xk)o - Ge), 


: {Ls (x) *r i 
* 
e [=x = (x) ,T 


s x= of “mS, on \* 
ae ety a) + (ant) She) ee: 
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# 
JT (x)dx = + Ty (x) ; 
duis) [rtGdx = 5 (O5@)-T2 0), 
* ul Tig) #1) 
JT (x) dx So dl aaa es ee 


interpretted as differentian formulae and used in recursion, give 


*! 





ale (x) 
th Aga h * * * 
Gh Eee 2(T, (x)4T on (xt. ; »+T, (x) ) +1 
(A-9) 
fe! 
1 Pan” 2(T. +I, Soin 
2 bot Bt ons GOT ea i te) 
Also 
ki Dele ne Oke 
(A-10) x T(x) = 1/2 A AGS ED) 
i=o ’ 
whence 
2k 
k 4 yy2k-1 2k, 
(A-11) x 5 1/2 2 Ce) 


* 
The orthogonal property of Chebyshev Polynomials is depicted by 


4 % | 1 (m=n=0) 
1” Tx (x) aT. (&) 
(A-12) f gx = ( 1/2 (m=n¥0) 
oO ¥x(1-x) 
0 (m#n) , 
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If a function f(x) is expanded as 
7 * * * 
ant == 
(A-13) f (x) 9 aot, (x) + aT, (x) ss ON ads aT) ; 


to evaluate it for any x in the range, the use of recurrence 


relation (A-6) gives > 

ay ea cS Gy ad 
2 ep 2 

where 


boo" a 4,4 bd a | + 2(2x-l)a_ 5 
(A-15) 


bo aEay + 2(2x-1)b 5 - b s(¢ = n-2(-1)0) - 


rrZ 


For special values of x, the results are simple and direct. 


id 
£(1) = 7 a + ay + a, + .t+ ay : 
uh 
(A-16) £(0) = 7a, 7 ay + a, ; 
Erne ee 
£(5) =7 4 a, ae ay, 
Also integration gives 
: 2 b 5 b tT" ts Pe) T" 
(A-17) IP SES Te ane ey Se og ae Eel 


oO 


—- ——) ee > 7 ial ae ee etd a, 
5 a | : , 7 
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where 
a -a 
r-l rt1 
be ie er Sel Lea e, 
eae) teh pte) , 
dvi ih. lode ERigy + (-1)"b 
20 ih 2 &} oy =) n+l 
Thus 
1 
(A=19) [> £(x)dx = 2(b,+b,1b,+...) 
fe) 
ee nt ie bad Ug 
20 bs3h G2 3054 yey) ia ras 
For differentiaition of f(x), if 
(A-20) £°(x), = 1a ee a. le +a! ; ’ 
2 el: n-l'n-1 
and 
(A-21) fs) — a"T" + aT" + + gt! * 
2 00 lige: : n-2° n-2 - 


where primes denote differention w.r.t. x, we have 
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so that ultimately 
_ ' 
(A-23) a. 4[(rtl)aj+(rt3)a ate. ’ 


so far as the terms vanish. Similarly, for the coefficients of 


second derivative 
(A-24) ay = 16[(r+1) (rt+2)a_ jt2(r+2) (r+4)a 
+ 3(r+3) (rt+6)a ot...) 


so far as the terms vanish. 
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APPENDIX B 


LISTINGS OF FORTRAN IV AND APL PROGRAMS 


The following pages consist of source listings of the 
FORTRAN IV programs GAS1, GAS2, subroutine CHSMY and APL 
program GAS3, subprograms LTINV1, LTINV2 and CHSMPY. The 
subroutine CHSMPY multiplies two Chebyshev Polynomials. The 
subprogram LTINV1 computes Ren from ae ; subprogram 


ttl 
LTINV2 computes 2 from ae Before calling GAS3, the values 


rel 
of A(=a,), Al(=a,), A2(=a,), TI (= initial temperature TW» Dy 
RHO(=p ), CRHO(= C)s K(=k), LSMAL(#2), MSMAL(=m), DELTA(=A), 
OBSVS(= time steps) and DTHETA(=A8) must be specified as global 
variables. V specifies the vector of two elements, Ax and 
distance-steps respectively. The output in each case are 


matrices Co: Cy +e M x Cc, and Ty where time increases along 


rows and distance increases along columns. 
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ITA bas YM2HD sotrucidue ,.S2A0 ,L2A0 emaxgorq VI WARTAOY 
edT .YOM2HO bas SVMITI ,IVATTI emesgoxqdye .€2AD mergozq 
sat .eleimonylol vereydsdd ows estiqtiium YIM@H saksvozdus : 


" Mexgoraqde ; “A moxz aA _esjuqmos IVMITI mexgozqdue 


i ssiuqmoo SVWMITI ' 


a (7 ezusjsisqmes LIetiint *) iT « (o8=)SA «(sea Cana, to 


“ geukev edd .€2AD gotileo sxotea ae moxt 


e(Se)ATING , (m=) IAMOM .(2=) TAM2I . (of=)a C9 =)OHFD .( 9=)OHT 
isdolg ae batitosqe sd tevm (94=)ATSHTG bos (eqete smts =)evedo 
bes xh ,stasmsie owt to toJosv siz esittosqe V .esidsizev 
938 sas> ross at tuqivo siT .ylevisoeqest eqoda~sonsteth = 

gools sevxestont amis sxsw oF bas P x M $9 tg esolztsm ‘ vz 


-aamylos gnols asessxont gonsyetb bas ewor — y 





. ss 


PAPA OI OS 


aD 


GAS) 
GAS ABSORPTION WITH CHEMICAL REACTION 
COMPUTATION OF ALL COEFFICIENTS RBIN) 
PREVIOUS, GOEEEICTENTS PRIN) TC. 


INOPUT DATA IS SCALED BEFORE INPUT 


BOI IIR RII ek a gtk ik aka kok kote ak ak kak ake ak ake ake ofc afe ate ak 
*BOUNDRY CONDITIONS 

% FR=1.0-71 .S¥*THETA 

* GR=2464789%*THETA 


ARIK ORR ORR IO Fok tek ica ak oak ak ak ak ge aka ok ak ak ake eae 


UGUBLE PRECISION RB(100)>PB(100) 
RG(200),PG(200),RP(200) ,PP(200 
RU(200),PU(200) ,R8V(200) ,PV (200 
FR{100),GR(100);F(100) ,G(100), 

25),TEMP0(25,25), ALPHA,BETA,GAM 
At sA2,9Dy,DELTA, KC, RHO,CRHD,TOLC 
RLyRM,PHI,PHIBAR,PEROR 

LASTR=4 

ISTEP=10 

X=4.50~-7 

DTHETA=1 .0D~-6 

READ(5,30) AyAisA2,sD,DELTAsKC,CRH 

READ( 5230) TO,LC, MC 


OW Wo ND et 


FORMAT (3010.3) 
X INC =X 
NN=1 


BETA=1/(D*DTHETA) 

ALPHA=BETA+LC/D 
GAMMA=RHO*CRHO/ (KC *DTHETA ) 
LAMBDA=DELTA*LC/KC 
RFOOTH=ALPHA**0.5 
NU=GAMMA**Q .5+ROOT 

WRITE (6,33 ALPHA sBETA,GAMMA,LAMB 
FORMAT (' LALPHA= Lee OA 
"CAMED A=", FT 2. 5s OX, *NUR*, E12. 5) 

Wee he Won oe) He rae LAS TAs Xie lo heP 
FORMAT (CY DIHETA=" ,FLO.9, ! TIME=-S 
i oe 6 eae BDISTANCE=STEPS=" 412) 
WRITE (6,236) AG AL Rees DEL faaKhG, Ge 
FORMAT(!? A="',E14.6,' AL=",E10.3;,° 
ey OE LTAS',F610.3 5!) KG="5F10.3547GRH 
WRUGLE(G, 327) TOsLCy MG 

FORMAT (*OINITIAL~TEMPERATURE=', E} 
os) 


bet 


COMPUTE CONCENTRATION AND TEMPERATURE 


ETC. 


Ae 


KA 


rXH(200),PH(200), 
),8Q(200),PQ(200), 
dyRW(200),PW(200), 
CONO(25,25),CONO1 (25, 
MA sLAMBDAsROOT,NUSAy 
Oy Leable Me sDTHETAY Xs 


V~RHO 


DA, NU 
SP ele. Fs OX y GAMMA] 47,622.55 ,555 
LEP Sy lee OLS TANGE-SiIL PS bie 
HO,RKRHO 

AZ EO. D=' sELO.S,y 


B=" pe toe sy RHOS4 5,610.39 


Oeds’ LSMAL=", 610.2," MSMALS*, 


AT THE sSURPACE 
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ee : 4 : a , 7 2 1h) ai Ty fkS 
| Peer i. pT (l | rsanely Paete 
ty iio eee 
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Poke 
iv. (ret are 7 
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ao i | 
Irina = 
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LATS 1 Gee we eas 7 
AV) JOee Ls oA AS i 
Ue ee peer irre 
ridge ee, OF Se ae 
UM s SORA Jy At) Ate A Donne ae ee 


PSU aB ity) WAMBAD! (RE ge aS 3d RTO TE rhe, Sete 9 ee 


wast Pe 


(2.513 4 Veli eke + PRL By = ORR oe 
~e oe 7 - : fi, Ge M4 i «> VTi ae 


L2 -Sa7e* @34AT Cin WELEteah Tyan - hi tet taira ath oe | 
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ea Peas Ane} Ge saan ler (21 Tet 
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40 


ee) 


PCa) eG ike) 


45 


150 


DO 40 1=1,LASTR 
GF CT )=TO#AX] *OTHETA 
FRCL P=AL+A2*GR (1) 
F(T)=FR(T) 
G(I)=GR(T) 
CONTINUE 

THE TA=OTHETA 

M= il 

R=0.0 

M2R=1 

M2R2=1 


old) Fl 0 
RM=0.0-LAMBDA*FR(M)/(ALPHA~GAMMA) 
RL=GR(M)-&M 

RP(i)=RB(1)*RG(1) 
RQ(1)=RB(2)*RH(1) 
RVCLI=ERL¥RP(1)/(NUSNU~ ALPHA) 
PWCL)=RM*RQ(1)/(3.0*ALPHA) 
RU(L)=O.0-RV(1)-RW(1) 
PH(L)=(D*PIXAL PHA) **0.5 
PHIBAR=0.0 


OUTPUT 
PA er el) eG Th 275 


TESTNEGR LASTR. SET PB=RB ETC. INCREMENT. 
CONTINUE 
TF(M.sEQ.LASTR) GO TO 300 
DO 150 I=1,M 
PB(I)=RB(1) 
PH(I)=RH(1) 
PG(I)=RG(T) 
PP(T)=RP(T) 
PO(I)=RQ(T) 
PU(I)=RU(T) 
CONTINUE 
DET SS: L=1 M2 
PV(T)=RV(T) 
PW(T)=RW(T) 

CONT TINUE 
THETA=THETA+DTHETA 
M=M+] 

P=M—] 

M1=M-] 

M2=M~2 

M2R=2%*M~] 
M2R1=2*M=2 
MZR2Z2=2*M+ 3 





G COMPUTE “RB 

FOCM)=BETA*PR(M=-1 )/ (2.0*R¥ROOT) 

PPAM.EG se)” GO TO f70 

DO 160 N=2,M] 

FN=N 

RB(NJ=RETA*PB(N-1)/(2.0*(EN-1.0)*ROOT)+ENFRBI(N4+1) /(2.0*RO0T) 
160 CONT I NUE 
170 RBCL)=FR(M) 


C 
c COMPUTE. RH 
RH(M)=1.0 
RH(Mi)=2.0%*(1 .0-GAMMA/ BETA) *R*ROOT/ (ALPHA-GAMMA) FRB (ML) /RB(M) 
TE(M.EO.27"°°GO TC iso 
DQ 180 N=1,M2 
EN=} 
RH(N)=0.0-(EN+] .0) ¥EN*RH(N+2)/ ( ALPHA~GAMMA ) 
j +2 .0*ROOT*EN*RH(N+1)/(ALPHA~GAMMA) 
2 —240*GAMMA*ROOT#R*XPH(N) /(BETA*( ALPHA-GAMMA ) )4+RB(N)/RB(M) 
180 CONTINUE 
C 
G COMPUTE. RG 


190 RG(M)=1.0 
[PUMeEQ.2) 50 TO 210 
DO 200 N=2,Ml 
EN=N 
RG(NJ=R¥*EPG(IN=-]2 J/CEN~1-0)4EN*RG(N41)/(2.0*GAMMA **0.5) 
200 CONTINUE 
210 RM=0.0-~LAMBDA*RB(M)/(ALPHA-GAMMA) 
RL=GAMMA**0.5*RL/(2.0%R) 
RG(L)=(GRI(M)-RM*RH(1))/7RL 


C COMPUTE RP,RQ 
DO 230 N=1,M2R 
RP (N)=0 
KQ(N)=0 
DO 220 J=1,™ 
Peele ait de Ou 1g 2 
Peo oU Ltr 220 
RPIN)= RP(N)+KB(J)*RG(IN41~-J) 
RQ(N) =RQ(UN)+RB(J) *RH(Nt1~-J) 
220 CONTINUE 
£25 FORMAT(20X,' N=" ,110.5X,*P(N) =',612.5,5X,* QIN) =',E12.5) 
230 CONTINUE 


c COMPUTE RV,RW 
RV(M2RY=RL¥KP(M2R)/ (NU*NU-ALPHA ) 
RW (M2R)=RMERQ(MZK)/(32.0* ALPHA) 
RV(M2R1)=4.0*NU*K*ERV(M2ZP) /(NU*XNU-ALPHA) 
j #RL*RP(MZRI)/ (NU*NU~AL PHA) 
RW(M2R1)=8.0*ROOT#R¥ERW(M2R)/ 03.0% ALPHA ) 
i tRM*ERQ(MZ2R1)/(3.0*AL PHA) 
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DO 240 N=L yaMZ2R2 
FN=N 
RV(N)=0.0-(CEN+1 60) *EN®RV(Nt+2) /(NURNU~ ALPHA ) 
+2 .O*NU*E NERV (N41) /(NUSNU-AL PHA) -BETA*PV(N) Z(NUENU -ALPHA) 
tRLERP ON) J/(NU*MNU- AL PHA ) . 
RWON)=0.0-(EN+1.0) ¥ENSRW(N4+2)/(3.0*%ALPHA) 
+4 OR IOTRENKRW(N41) /(2.0%AL PHA) -BETA*PWIN) /(3.0%ALPHA) 
+KM¥RQ(N) /(2.0* ALPHA) 
CUONTLINUE 


COMPUTE RU 


RUC(M)=BETA*PU(M-1)/(2.0*R*ROOT) 

Pai WeG eee eGU 10 260 

DO 250 N=2,M1 

EN=N 

RUCN)=BETA*PU(N-1)/(2.0*(EN-1.0)*ROOT )F+EN*RU(N41)/(2.0*ROOT) 
CONTINUE 

RU(1)=0.0-RV(1 )-RW(1) 


COOMEUte ePitl.s..~PHIBAR 


4 
- 


t— 


PHI=((D¥*¥PI*¥(R4+1.0))**0.5) *(ROOT-KRB(2)/FR(M)) 
PHIBAR=((D¥*PI*(R+1 .0))**0.5)*(ROOT*RW(1L)4RV(1) *GAMMA*#*0.5 
Sete =e Vt 2h Wl 7K (Mo) 

FACTOR=(D*PI*(R+1 .0))**0.5 

TERM=ROOT-RB(2)/FRI(M) 

FORMAT (5X _'***ee*XTRACE FACTOR » TERM = *,27F12.5) 

CONTINUE 

BSUM=0.0 

GSUM=0.90 

HSUM=0.0 

USUM=0.0 

VSUM=0.0 

WSUM=0.0 

BSUM=BSUMF#RB(T)*&X**(T-1) 

HSUM=HSUMF#RM *PH (1) *X** (1-1) 

GSUM=GSUM+#RL ERG( I )*X**( 1-1) 

USUM=USUM+RU (I )*X**(I-1) 

CONTINUE 

DQ 290 [=1,M2R 

VSUM=VSUM#RV (T )*X*#*( 1-1 

WSUM=WSUM#RWCT )*EX** (I~) 

CONT [NUE 

CO=BSUM*XDEXP(0.-0-ROOT*X ) 

TO=GSUM*DEXP(0.0-GAMMA**0.5*X) tHSUM*DEXP(020-ROOT *X) 

CL=USUM*D EXP (0.0-ROOT*X) #VSUM¥DEXP (0 .0-NU*X) 
+WSUM*DE XP(02.0-2*ROOT*X) 

CONO(NN»yM)=CO 

CONOL(NN,M)=CO+(MC/D) ¥Cl 

TEMPO(NN,M)=TO 

FR(M)=CO 


) 
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GR(M)=TO 


PPO MeEO. LJ GO TO 145 


GO TO, 100 
NN=NN+1 


LECNNGGE. ISTEP) «GQ TOn40 


PEROR=1 .002* ( 
CONO( ISTEP 
WRITE(6,500) 


CONOM(ISTERSLASTRI=CONOLISTEP,LASTR))/ 
PEROR 


FORMAT(*O PERCENTAGE-ERROR= 4,025.15) 


WRITE(69127) 
FORMAT (# | kx % 
WEITE( 6,226) 
RUEMAT (4s «SURE 


CQ eH Y 
(FCT ),1T=1,LASTR) 
feet (1X oA tek 5)a) 


DG S310, Lalo: bSTEP 


WePrE CG Pe) 
FORMAT(*' M=*, 
WRITE(65128) 
FORMAT ("1 x kx 
WRITE(6,126) 


v(CONO(T, J) pJ=LyLASTR) 
Ais StleX » D250 1 Si) 


CO+M*¥C] MA HY 
(FAT bie Teles tol Spey) 


DO 9340) sleds LSE? 


WRITE(6,123) 
WRITE(6,124) 


Let CONOL Olio: Jae hah pLAS TRG 


FORMAT ( © 1 Ak ae otc oe te 2 ok Xe ae a a ok ake ok ake ‘Bie 3 RC RO Rok ok doko deka ak ack koakogeok ack @) 


WRITE(6,126) 


(G(T ),1=1,LASTR) 


DO, 320) d= 14 bSTEP 


WRITE(6,123) 
ST OP 
END 


Liowt We MEO Cl he) sa, LAST R.) 
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GAS2 
PROGRAM TO CALCULATE THE CONCENTRATION OF GAS AT ANY DEPTH 


CONSTANTS USED ARE AS FOLLOWS * * * * % 

PCL Teeete MenoOURCI IN Bel.Ue PER the MOLE J 

RH Ss" MOLAR DENSITY 

Cine. SULA seco l Pic HEAT. 

Kee rteeMmaAtl DENSITY 

(2e ees es felt iae” ToMPERATURE “AT THE SURFACE 
Weeteose ore USl VI TY 

tities tO tae TIME OF OBSERVATION 

LSMALS REACTION RATE IS K(T)=L4#M.T WHERE L=LSMAL 

MSMAL =: AND M=MSMAL. 


INTEGER OBSVS,0OBS4,0BRS1,0BS21 

DOUBLE PRECISION D,ALFAL,KyDELTA,RyT IME, LSMALy, MSMAL oL 
] MyNU,yLAMDA CO, TO,C1LsASUM,RSUM,CSUM,USUMy VSUM,WSUM, 
2 SALFA,SGAMA,NU,P,PEROR,AR(100),¢RR(100),CR(100), 
3 AP 1(100),8R1(100),CR1(100),UR(100),VR(200),WR(Z00), 


4 


UR1T(100),VR1 (200) WRI (200),LR1i( 200) »MRL(200),F(100), 
5 G(100),FF(100),6G(100) ,CONO( 25,25) ,TEMPO( 25525), 


6 CONO01 (25,25) ,TOs,sA,Al-A2-RHO;CRHO 
DOUBEE PRECISION CONT( 25,25) 

FORMAT (8010.3) 

READ(5 +1 )A,A1+A2+D,DELTA,K,CRHO,RHO 
READ(5,1)TO,LSMAL »MSMAL 


Per PAPI SeE THE * CONSTANTS 
Seber VESVS-TIME-STEPS AND [STEP=DESTANCE=STEPS 
Be miuey P=-Dl STANCE STEP=SIZE, DIHET =TIME STEP=STIZE 
NBSVS=4 
TSTEP=10 
Orie (a1 emt 
R=4.50-7 
ALFA1=K/ (RHO*C RHO) 
L=LSMAL/D 
M=MSMAL/D 
BETA=1/(D*DTHET ) 
ALFA=BETA+#L 
GAMA=1/(ALFAL*DTHET) 
LAMDA=DELTA*LSMAL/K 
SALFA=SQRT(ALFA) 
SGAMA=SQ2T(GAMA) 
DSFG=SALFA-SGAMA 
NU=SALFA+SGAMA 
WRITEC6 2) ALFA,BETA,GAMAsLAMDA, NU 
FORMAT (* ALFA= Sy Re Ale ee BETA= Pe decreOren. GAMA= vy, 
i ee ose CAMDA=";DIi2~+6;' NU=",DL2.5) 
WRITE(6,3) DTHET »,OBSVS,R,ISTEP 


mAs 


Uy i 
wd 
’ Sys 
Sere 
’ Ne Ai TA y 
; ; 
} 1 
J i 
= i 
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3 BORMAT(’ DELTA-THETA= *,D12.7,' TIME=STEPS= °,15 
Ul MiPoeeTA=X%=",Dl2.7)") DISTANCE=STEPS=! ,15) 
WRITE( 46,4) Bly A stAes Dy DELTA SK 
a FORMAT (! Nate Ol Sy 5 ’ VIR Suey aye an Ag sips. .5 7 D=! 
l yORo com MIDE PA= ODL Anh) ahi yO2S 15) 
WRITE(6,5) CRHO,RHO,LSMAL »MSMAL ».TO 
5 moments CRHUS'», 015.5 ,% Ph eg Lt ie vy LSMALH=!', 015.5 
ee es Somat", 015.5)" TO=",015.5) 
OBS4=OBSVS+4 
NO 10 N=1,0BS4 
AK(N)=0.0 
BRIN)=0.0 
CR(N)=0.0 
UR(N)=0.0 
ARIT(N)=0.0 
ARL(N)=0.0 
CRI(N)=0.0 
UR?(N)=0.0 


9 


10 CONTINUE 
OBS21=2*0BSVS+3 
DO Pe ON =-1 5 OBS21 
VR(N)=0.0 
WRIN)=0.0 E 


VR1L(N)=0.0 
WRL(N)=0.0 
LR1(N)=0.0 
MR1(N)=0.0 


me CONTINUE 

CPs ArmeGAUGUMATIONS FOR ONE STEP IN "TIME? DIRECTION 
NN=1 
D0 12 T=1,90BSVS 


G( 1 )=TO+A*I*DTHET 
PANY =AL+A2*G(T) 
PRQL)=FC13 
GG(1)=G(1T) 

ed CONTINUE 
AR(1)=2.90 *F (1) 


CR(L) = —2*LAMDA*F(1)/ (ALFA - GAMA) 
BRIG H=220FOGHIO-ERIT) 

VR (1) = F(1L)*BR(1L)/(NU*NU-ALFA) 
WR(1) = F(1L)*CR(1)/(3.0*#ALFA) 


UR(1L)=-(VR(L)+WRO1)) 

CO=AR(1)/2-e0*DEXP(0.0-R*SALFA) 

TO=BR(1)/2 e0*DEXP (0 20-R¥*SGAMA )#+C2(1)/2.0*DEXP(020-R¥*NU) 
C1=UR(1)/2+ OXDEXP(0.0-R¥S ALFA) #VR(1)/2.0*DEXP(0.0-R*NU) 
1 +WR(1)/2.0 *DEXP (0.02. 0*R¥*SALFA) 

CONO(NN,1 )=CO 

TEMPO(NN,1)=T0 

CONOL(NN,1 )=COt(MSMAL/D) *C 1 

F(10=C0 

G(1)=T0O 
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rap 
*4=02 









1=2 
20 CONT LNUF 
no Pus elee Bane J 
N=I+2—J 
ARLINDSARLIN#2) + (2*N/(RESALFA))FARL(NGL) + (R¥BETA/ 
1 (8*SAL FAY) *( (CAR (N=21)-AR (NFL) /Z(N=1)=(AR (NFL) -AR (N43) /(NF1L)) 
CONTINUES 
AR1(1)=0.0D0 
pa 24 N= 2,1 
P=(-1.0000)**N 
ARL(L)=AR1(1)+ARi (N) *P 
24 CONTINUE 
ARL(1)=2.0 *(F(E)#AR1(1)) 


J 
NM 


DO 2Ba iJ 2h344 
DUR eS ee ae eT es: 
N=I+2~-J 
BRI (MN) =BRL(CN#+2)+2*N¥*BRI(NtL)/(R*SGAMA )+4R¥*¥SGAMA¥((RR(N-1) 
€ ~BR(N4+1))/(N-1)—(BRI(N4+1)-BR (N43) )/(N41))/39.0 
28 CONTINUE 
DOSS 00JeA1 7 I 
N=I+1-J 
CRI(N)=8*N*SALFA/ (R¥( ALFA-GAMA) )*®(CRI(N4t1)-CRL(N4+3))+ 
1 CRLINt2)*(14+16*N* (NFL )/(R¥*¥R*®(CALFA-GAMA) ))+(0N/ (N42) ) * 
2 (CRLIN#2)-CRLIN+4) )-(GAMA/ (ALFA-GAMA ) )*((CRINI-CR 
es (N42))—-(N/(N4#2))*(CR(N#2) -CR (N44) ))-(LAMDA/CALFA~GAMA)) 
“& *( (ARLIN) “ARLIN+2) )—(N/ (N42) )®(ARL(N#2)-AR1 (NF4))) 
30 CONTINUE 
BR1(1)=0.0 
NQ 32 N=2,I 
P=(-1.0000)**N 
Bed CMe aBR1 C1) + (BRIUN)4+CRi(N))*P 


a2 CONTINUE 
BRL(1)=2.0 mUG(1)4+8R1(1})-CRI(2) 
(eo Atiie CHSMY Ci, Anis lis BRLs CRY) 
Gta os IGHSAY (I1,ARL,I1,CR1,MR1) 
T2i=2*I1-1 
DN 38 J=2,!I 
N=1I+2-J 
URL(N) =2*N*URL(NtL)/(R*SALFA) + URL(N42) + (BETA*®R/( 8% 
L SALFA)) *( (UR(N=1)-UR(Nt1))/(N-1)- (UR (N#+1)-URON43))/ 
2 (N+1)) 
38 CONTINUE 
DN 40 J=b, 121 
N=2*I~-J 


VRLINJEVRIUNE2 (1-1 GNF (N41) / (RER*®SGAMA% (NUFS ALFA) ) ) + 
( B*ENXNU/ (R*¥SGAMA*(NU+SALFA)) Y*¥(VRLONtL)-VRLONt3) + 
(N/(N4#2)) ¥(VR21ONF2)-VRLINt4) = (BETAS (SGAMA® (NUFSALFA 
Y)*CCVRONIAVRON+2) 2) = ONZ(N#2)) #(VRONF2 2° VRONF4)) D4 
COULRLON)J-LRI(Nt2D)-(N/(N42)) *(LRIC(N+2)-LRIUNt4S)))/ 


(SGAMA*(NU+SALFA)) 
WRL(N) =WRLCN#2)*( 1-1 OEN¥X(N#L) /(3*R*R*ALFA)) + (16*N/ 
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40 


44 


48 


ae. 


t~s 


Ww 
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(34R*XSALFA))*(WRI(NFL)—-WRLIN4+3) )FON/ (N42) ¥(WRIINH2 
J-WRLINE4 ) )- (BETAS ( B#AL FA) ) ¥( (WR ON) =WRONt2) - CNS ON42 
))#(WRINt+2)—WRONtG)) DFC (OMREON) @MRLCNt2))-(N/ (N42) 7 * 
(MR LINF2 )-MRLIN+4)))/03*AL FA) 

CONTINUE 

WEL es ? NE al 

P=(-1.00D00)**N 

ee Cr) aur CL yu Rt UN JP 

CONT INUE 

DO 44 N=2,I21 

P=(-1 .0D00)**N 

VWRLAEVWRLECVRIUCN) tWRIOCN) ) ¥P 

CONTINUE 

UR1(1)=2.0 *(UR1L(1) + VWR1) 

N90 46 N=1,I 

AR{N)=AR1L(N) 

BR(N)=BRI(N) 

CR(N)=CRIIN) 

URIN)J=URLIN) 

CONTINUE 

DO 48 N=lyt2 

VR(N)=VRICN) 

WR(N)=WRIOCN) 

CONTINUE 

ASUM=AR1(1)/2.0 

BSUM=BR1(1)/2.0 

CSUM=GRi(1)/250 

USUM=UR1 (3) 72.0 

VSUM=VRL(1)/2-0 

WSUM=WR1(1)/2-0 

Do 50 J=2,1 

ASUM=ASUM+AR (J) 

BSUM=BSU MtBR (J) 

CSUM=CSUM+#CR (J) 

WSUM=USUM+FUR (J) 

CONTINUE 

COH=ASUM*DEXP (-R*¥SALFA) 

TO=RSUMEDEXP (-R¥SGAMA) +C SUMEDEXP( ~R¥SALFA) 

DG. 52sN=2 > 121 

VSUM=VSUM4VR1 (J) 

WSUM=WSUM4tWR1I (J) 

CONTINUE 

C1 =USUM#DEXP (-R&SALEA) #VSUM*DE XP (-R*NU) WSUM#DE XP (~2 0% 
R*SALFA) 

CONO(NN,1I)=CO 

TEMPO(NNy,1)=TO 

CONOL (NN, 1 )=CO+ (MSMAL/D)*C1 

FE (1)=CONO(NNyT) 

G(1)=TO 

[=I +1 

Tele wis y oJ 00 Zz 


{ { :? ‘gem a9 


Mo EO 


Ws FO) "4 
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L100 CONT TINUE 
NM=NN41 
PEONIVSLERISTEP) I GDITO 12 
PEROMR=1.0D02* (MSMAL/D)*CONLCISTEP,OBRSVS)/CONO(ISTEP,OBSVS) 
WRITE(6,6) PEROR 


é FORMAT("O PERCENTAGE-ERROR= ',0D25.35) 
WRITE(6,7) 
7 FORMAT (Nir CQ okxn) 


WRITE(6,15) (FFECI),1=1,0BSVS) 
09 54 I=1L,ISTEP 


b4 Me lyeC6,3)71L9(CONO(1 ,J) ,J=1,08SVS) 

rs} RORMAP{® XSTEP=",125401X%,D25.15)) 
IRITE(G)9) 

9 FORMAT( PY) kee 860 ck KH) 
Hees, 15) (6G(1),1=1,0BSVS) 

15 mennart' SURFACE ',4(1X%,025.15)4 
DO S6° f=? STEP 

56 BElVEbGs S9eh 1, VIEMPOLIL.J),J=1,0BSV5S) 
WRITE(6,16) ; 

16 FORMAT (*1***e COF+(MSMAL/D)*¥C1L «ex! ) 
POIG0el=bsI STEP 

60 Wet e 665108 4) DOG UNOL (1,3) ,.J=1,CB8SV3) 

S EGP 


END 
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(a) 
™ 


THIS SUBROUTINE MULTIPLIES TWO SHIFTED CHEBYSHEV POLYNOMIALS 


SUBROUTINE CHSMY (MyAyN,BeW) 


ChOgee (Ota ee) AN) (BO ett Ble liiteece! 
YIELD THE PRODUCT (WO/2.TOTtWIL 11+...) OF 
Cee KK TNPUTEVECTORS (AOsALy eee) & (BOrBly eee) 
C***KXXOUTPUTHEVECTOR (WO,Wly eee) 
C 


iat 


REAL*8 A(1i0075 8 (100), 01200) 
FAR (ER ey 9 Gat lead By Ar ae 
B{1)=B8{(1)/2.0 
K=M+N-] 
ope LOe Jes K 
Ww{J)=0 
CONTINUE 
Ni=l1 
AA(J)=A{(J)*BCN1) 
CONTINUE 
NO 14 [=1,M 
INL=T4#Ni-l 
WC INL) =WCUINL )+AAC(T) 
CONT INUE 
(ee (MAGT.NL? GO" Tb 20 
Wad ey J=1,M 
I=N1-J+tl 
W(T J=W{ 1) +RAL J 
CONTINUE 
NL=N1 +1 
Peni Nove hy GU TO 44 
GO. TOf 30 
Do 22) J=1 NL 
[=Ni-J+1 
W(L)=W(1)+AA( J) 
CONTINUE 
M1=M-N1 
DO 24 J=1,M1 
J1l=J+1 
MI t=N1+J 
wW(J1)=WCJL)F+AA(N11) 
CONTINUE 
Ni=N1 41 
ire ONYOLER NOE GOD) it 
CONTINUE 
00 322 IT=29K 
Wil J=W 1) /2 20 
CONTINUE 
RETURN 
END 
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Vecee omy. CONDO CONIA TEMPO TMARSIMER*IMCR: TIMERS IMER AR: 

ARLODRABRAs On SCR SURTURLI URSA. VA SVAAVEASWO SWISS WRAs ES 
Be Gua 2. fel 2 

C1) FF«GG+F*G+OBSVSp0xI+J+N+l 

D2 GAO Sl sCGGelI«GlLTIeTIOrAxTxDTHETA 

ee RELL NGE 1 |tA ox Ax xO THE A 

C4] >((T+I+1)SOBSVS)/GAS31 

Ss | GC Aa eA ats Worx KA ee es SS Kets Ket DEL 
TAS). DELIA ) 

6) tot 


re CYRHO="*RHO:'* CRHO='3CRHO;'° LEMAE=" 3 05MAL; '* MOMAL=",; 
MoMA st TO=! : 7 TO) 
[8] if. Q 


bo) CLALPEA 
E20) (GLBETA 
C11] ('GAMMA 


'sALFA<«(LSMAL+1+DTHETA) =D) 

Vanna se OxD ee AD) 

'.s GAMA*RHOxCRHO+(KxDTHETA) ) 

£12] ('LAMBDA '; LAMDA*+DELTAXLSMAL+K) 

bret ACI DELTA-TAETA '; DTHETA) 

C14] ('TIME-STEPS '; OBSVS) 

eo) RUDE LTA RX '.DX+V[1]) 

Pipe eee me ANCE =o 0k bo: rye |S) 

Si ae LOLAG=T IME ', OBSVSxDTHETA) 

[18] ('fOTAL-DISTANCE Per alc VL 21) 

19 ! ' 

a CON0+CON1+TEMPO<+(V[2],0BSVS)p0 

[21] IMAR+(1 1)912¢EPS1«DXx(SALFA*ALFA* 
Ge aay, 

[22] IMBR<(1 1)912EPS2<«DXx(SGAMA*+GAMA*x 
Os Ds 2 

Poag) 2MCR+(2 21 2X11) ,0,(-ePst:(xATi*2)),(1tXl1<( GAMA- 
ALFA)x(DX*2) #16) 

C24] IMER+*LTINV2(2 yo (14X12) 70, C-CEPS Tree s2) tts i 2s 
2)),1¢XI2<-(( GAMA+2xSALFAXSGAMA ) x(DX*2)+16) 

Fol FMFR<IPINV2(2 2)0(1tX13) ,0, 0-28 PS12(XT3*2))  (itArse— 
3xALFAx(DX*2)416) 

[26] GAS32:IMAR«LTINV1 IMAR 

[27] IMBR+LTINV1 IMBR 

p38] IMCR*LTINV2 IMCR 

29] IMER<LTINV2 LTINV2 IMER 

[30] IMFR+LTINV2 LTINV2 IMFR 

£31] 3((J+J+1)<OBSVS)/GAS32 

Pee cenGAo ou std el 

fea) ARs, 2@F Ul I 

[34] CR+,-LAMDAxAR+(ALFA-GAMA) 

Pasa vesre,2<Chl) =ck 

[36] VR+«,BRxFLIJ2( (NU*SALPA+SGAMA )*2)-ALFA 

(37] WR+e,FLIJxCR+3 

[38] UR+<«,-VR+WR ; 

[39] CAS34:CONOLN s7I*+FPLII<((+/AR) -ARD1]22)x*-DXxXSALPA 

[40] PEMPOEN sf GL I<( ((+/BR) -BRL1122) x*-DXxSGAMA)+((+/CR) 
~-CR(1]22)x*-DXxSALFA 

Cae) +CAS35 , (URX*MSMALXUR#D) , (VRX*MSMALXVR#D) , (WRXMSMALX 
WR+D) 
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